function [Y,Xx,Ea,Eu,P] = mysimulateper(This,S)
% mysimulateper  [Not a public function] Simulate period by period.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

%--------------------------------------------------------------------------

ny = size(S.Z,1);
nx = size(S.T,1);
nb = size(S.T,2);
nf = nx - nb;
ne = size(S.Ea,1);
nPer = S.NPer;
nAlt = size(This.Assign,3);
iAlt = min(S.ILoop,nAlt);

% Main loop
%-----------
% Set time for nonlinear equation evaluation. We do period by period, hence
% the period being evaluated is always tt=2 (after init cond included).
tt = 2;

% Preallocate output data.
Y = nan(ny,nPer);
Xx = nan(nx,nPer);
Ea = S.Ea;
Eu = S.Eu;

if S.IsUpdate
    eqtnN = S.Update.EqtnN;
    posUpdate = S.Update.PosUpdate;
    nPosUpdate = length(posUpdate);
    P = nan(nPosUpdate,nPer);
end

testFn = @(x) ~any(strcmp(x,{'off','none'}));
isDisplay = ...
    (S.IsUpdate && testFn(S.Update.sstate.optimset.Display)) ...
    || ...
    (isequal(S.Method,'global') && testFn(S.OptimSet.Display));

if isDisplay
    nDec = floor(log10(nPer)) + 1;
    textfun.loosespace( );
end

for t = 1 : nPer
    if isDisplay
        doRptSolving( );
    end
    
    S.y = [ ];
    S.w = [ ];
    S.x = [ ];
    
    % First-order approx simulation first.
    S.Alp0 = S.U \ S.x0;
    
    ixPer = false(1,nPer);
    ixPer(t) = true;
    
    S.M = simulate.linear.multipliers(S);
    switch S.Method
        case 'firstorder'
            [S.y,S.x,S.Ea,S.Eu] = simulate.linear.run(S,1);
        case 'global'
            [S.y,S.x,S.Ea,S.Eu] = simulate.global.run(S,ixPer);
        otherwise
            utils.error('model:mysimulateper', ...
                ['Function simulate( ) with sstateUpdate=true ', ...
                'not implemented yet for method=%s.'], ...
                S.Method);
    end
    
    Y(:,t) = S.y(:,1);
    Xx(:,t) = S.x(:,1);
    Ea(:,t) = S.Ea(:,1);
    Eu(:,t) = S.Eu(:,1);
    
    % Store current parameters and update sstate and solution.
    if S.IsUpdate
        P(:,t) = This.Assign(1,posUpdate,iAlt).';
        if t < nPer
            doUpdate( );
        end
    end
    
    % Prepare next period.
    S.x0 = S.x(nf+1:end,1);
    S.Ea(:,1) = [ ];
    S.Eu(:,1) = [ ];
    S.G(:,1) = [ ];
    if ~isempty(S.Anch)
        S.Anch(:,1) = [ ];
        S.Wght(:,1) = [ ];
        S.Tune(:,1) = [ ];
    end
end

return

    
    
    function doUpdate( )
        % doUpdateParams  Evaluate !sstate_update and update param values.
        if isDisplay
            doRptUpdParam( );
        end
        
        xx = [ [nan(nf,1);S.x0], S.x(:,1) ];
        xx(S.IxXLog,:) = real(exp(xx(S.IxXLog,:)));
        ee = [ nan(ne,1), S.Ea(:,1)+S.Eu(:,1) ];
        yy = [ nan(ny,1), S.y(:,1) ];
        pp = real(S.Assign(1,S.NameType==4));
        LL = S.L;
        % Absolute time within the entire nonlinear simulation for sstate
        % references.
        ttAbs = -S.MinT + t;
        
        % Evaluate sstate update equations
        %----------------------------------
        try
            D = eqtnN(yy,xx,ee,pp,tt,LL,ttAbs);
        catch Err
            utils.error('simulate:mysimulateper', ...
                ['Error evaluating sstate update equations at t=%g.\n ', ...
                '\tUncle says: %s'], ...
                t,Err.message);
        end
        
        % Update
        %--------
        D = D(:).';
        ix = This.Assign(1,posUpdate,iAlt) ~= D;
        if any(ix)
            if isDisplay
                name = This.name(posUpdate);
                doRptChange(name(ix));
                doRptUpdModel( );
            end
            
            S.assignpos = posUpdate;
            S.stdcorrpos = nan(size(S.assignpos));
            % Update sstate and solution.
            isError = true;
            This = myupdatemodel(This,D,S,iAlt,S.Update,isError);
            % Update Assign and solution matrices in S.
            S = myprepsimulate2(This,S,iAlt);
        else
            if isDisplay
                doRptNoChange( );
            end
        end
    end % doUpdate( )

    
    
    function doRptSolving( )
        fprintf('=====Period %*g/%*g. Solving equations.\n',nDec,t,nDec,nPer);
    end % doRptSolving( )

    
    
    function doRptUpdParam( )
        fprintf('=====Period %*g/%*g. Updating parameters.', ...
            nDec,t,nDec,nPer);
    end % doRptUpdParam( )

    
    
    function doRptChange(List)
        fprintf(' Change in');
        fprintf(' %s',List{:});
        fprintf('.\n');
    end % doRptChange( )

    
    
    function doRptNoChange( )
        fprintf(' No change.\n');
    end % doRptNoChange( )

    
    
    function doRptUpdModel( )
        fprintf('=====Period %*g/%*g. Updating steady state and solution.\n', ...
            nDec,t,nDec,nPer);
    end % doRptUpdModel( )
end % main
